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the  complex  modulus  with  the  use  of  fractional  calculus.  The 
use  of  fractional  derivatives  and  a  temperature  shift 
function  provided  a  complex  modulus  model  that  could  be  used 
in  a  finite  element  formulation  to  analyze  the  system 
response  of  a  structure  that  contained  both  elastic  and 
viscoelastic  members.  While  the  model  is  limited  to  low 
frequencies  and  use  of  the  complex  modulus  in  the  transition 
region,  the  system  response  is  a  relatively  simple  expression 
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Abstract 

The  purpose  of  this  study  was  to  demonstrate  the  use  of 
fractional  derivatives  to  capture  the  frequency  and 
temperature  dependency  of  viscoelastic  material  behavior.  To 
model  the  frequency  dependency  of  viscoelastic  material, 
fractional  derivatives  were  included  in  the  complex  modulus. 
Solution  techniques  were  performed  in  the  Laplace  domain  to 
allow  for  easy  manipulation  of  the  fractional  derivative 
terms.  To  incorporate  the  temperature  dependency  of 
viscoelastic  material  in  the  complex  modulus  model,  the 
method  of  reduced  variables  was  employed  with  the  use  of  the 
WLF  equation.  Wilh  the  frequency  and  temperature  dependency 
built  into  the  complex  modulus,  a  finite  element  formulation 
was  devised  that  incorporated  elastic  and  viscoelastic 
response  of  a  truss  structure.  The  formulation  was  limited 
to  the  use  of  the  complex  modulus  in  the  transition  region, 
the  region  where  the  damping  ability  of  viscoelastic  material 
is  maximized.  Quasi-static  motion  was  also  assumed,  which 
limited  the  response  to  low  frequencies.  0 

The  solution  technique  was  demonstrated  on  a  nine  degree 
of  freedom  truss  composed  of  aluminum  and  neoprene  rubber 
rods  subject  to  a  temperature  variation.  The  results  of  the 
example  problem  show  that  temperature  and  frequency 
dependency  of  viscoelastic  material  can  be  incorporated  with 
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the  use  of  fractional  derivatives.  The  simple  solution 
format  and  model  robustness,  along  with  the  theoretical 
basis,  provides  encouragement  for  additional  work  on  the 
complex  modulus  model. 


vi  i  i 


FRACTIONAL  CALCULUS  FORMULATION  OF  THE 


QUASI-STATIC  VISCOELASTIC  PROBLEM 


I .  Introduct i on 


1 . 1  Overview 

The  purpose  of  this  research  effort  is  to  demonstrate 
the  use  of  fractional  derivatives  to  capture  the  frequency 
and  temperature  dependency  of  viscoelastic  material  behavior. 
This  chapter  provides  a  basic  introduction  to  the  development 
of  the  use  of  fractional  derivatives.  In  chapter  two,  the 
fractional  derivative  is  incorporated  into  a  complex  modulus 
model.  The  effects  of  strain,  frequency,  and  temperature  are 
investigated  as  they  affect  the  viscoelastic  behavior.  With 
the  fractional  derivative  complex  modulus  model  established, 
a  two-dimensional  finite  element  formulation  is  developed  in 
chapter  three.  This  solution  technique  imposes  the  condition 
of  quasi-static  motion  and  incorporates  temperature  and 
frequency  dependence  in  the  transition  region.  The  use  of 
the  finite  element  formulation  is  demonstrated  in  chapter 
four.  This  example  problem  provides  a  solution  technique  for 
the  system  response  of  a  typical  space  truss  configuration 
composed  of  elastic  and  viscoelastic  members.  The  final 
chapter  is  devoted  to  conclusions  of  the  solution  technique. 
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1 . 2  Background 


Recently,  there  has  been  an  increasing  demand  to  model 
the  frequency-dependent  behavior  of  viscoelastic  materials. 
The  driving  force  behind  this  resurgence  of  viscoelastic 
behavior  models  has  been  the  increasing  use  of  viscoelastic 
material  in  engineering  structures  to  take  advantage  of  their 
energy  dissipation  characteristics.  One  common  approach  to 
modeling  the  viscoelastic  behavior  is  to  include  a  large 
number  of  derivative  terms  to  model  the  frequency-dependent 
stiffness  and  damping  characteristics  of  the  material 
(1:918).  The  result  is  a  complicated  model  that  produces 
even  more  complex  equations  of  motion  to  be  solved.  Another 
variation  is  to  store  the  real  and  imaginary  components  of 
the  complex  modulus  in  a  computer  (1:918).  When  the  data  for 
a  particular  frequency  is  required,  the  data  can  be  accessed. 
This  works  well  for  a  single  frequency,  but  this  process 
still  results  in  a  complicated  solution  for  transient 
analysis.  The  use  of  fractional  derivatives  to  model  the 
complex  modulus  has  shown  great  accuracy  and  simplicity. 

The  concept  of  using  fractional  calculus  is  not  new.  An 
article  by  Ross  (2)  provides  an  interesting  historical 
summary  of  the  development  ol  using  derivatives  of  fractional 
order.  Ross  states  that,  in  1695,  L'Hospital  posed  a 
question  to  Leibniz  of  the  use  of  a  derivative  of  fractional 
order.  Leibniz's  response  was  that  it  would  lead  to  a 
paradox,  but  he  also  concluded  that  "some  day  it  would  lead 
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to  useful  consequences."  Even  back  in  1819,  derivatives  of 
fractional  orders  were  discussed  in  a  text  by  Lacroix.  While 
the  idea  of  fractional  calculus  has  been  around  for  nearly 
three  hundred  years,  the  use  of  fractional  derivatives  seems 
non-existent  as  compared  to  integer  order  derivatives.  In  a 
series  of  papers  by  Bagley  and  Torvik  (1),  (3),  (A),  and  (5), 
the  development  of  the  use  of  fractional  calculus  to  describe 
viscoelastic  behavior  is  well  documented.  The  initial 
observations  by  Nutting,  and  the  subsequent  findings  by 
Gemant ,  Scott-Blair,  Caputo,  and  others  are  traced.  The 
fractional  order  derivatives  were  used  by  Bagley  and  Torvik 
(1),  (3),  (A),  and  (5),  to  model  the  frequency-dependent 
behavior  of  viscoelastic  materials.  The  advantage  of  the 
fractional  derivative  model  lies  in  its  simplicity.  By  being 
able  to  model  the  complex  modulus  with  as  few  as  three  or 
four  parameters,  a  very  accurate  least-squares  curve  fit  can 
be  obtained  to  model  measured  frequency-dependent  mechanical 
properties.  The  mathematics  involved  lends  itself  to  the  use 
of  either  Laplace  or  Fourier  transformations.  Results 
obtained  have  shown  accurate  descriptions  of  complex  modulus 
data  up  to  eight  decades  of  frequency  with  a  four  parameter 
model  (5:133). 

’rhe  use  of  the  fractional  derivative  in  the  finite 
element  formulation  for  determining  transient  response  is  in 
agreement  with  the  continuum  model  (1).  The  ability  to 
include  the  fractional  derivative  in  the  finite  element  model 
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provides  greater  flexibility  to  analyze  structures  containing 
both  elastic  and  viscoelastic  components.  In  a  more  recent 
paper  by  Bagley  and  Calico  (6),  the  fractional  derivative  was 
utilized  to  produce  equations  of  motion  capable  of  describing 
the  system  in  terms  of  its  initial  states.  The  construction 
of  the  state  equations  with  fractional  order  derivatives 
provide  additional  forms  of  feedback,  thus  enabling  engineers 
to  improve  feedback  control  systems. 

While  the  fractional  derivative  model  provides  an 
accurate  empirical  model  of  viscoelastic  material  behavior,  a 
theoretical  basis  would  increase  the  confidence  in  its  use. 
The  fact  that  the  fractional  derivative  formulation  created  a 
simple,  yet  quite  accurate,  representation  of  the  complex 
modulus  leads  one  to  question  its  fundamental  nature.  In  a 
paper  by  Bagley  and  Torvik  (4),  this  theoretical  basis  was 
established.  The  result  of  the  molecular  theory  by  Rouse 
"provides  us  a  non-empirical  basis  for  including  fractional 
derivatives."  (4:208)  Ferry,  Landel,  and  Williams  adapted 
Rouse's  theory  to  concentrated  polymer  solutions  and  polymer 
solids  with  no  crosslinking  (4:208).  The  result  was  that  the 
viscoelastic  shear  stress  was  a  function  of  the  fractional 
derivative  of  the  shear  strain  history. 

The  successful  ability  to  model  frequency-dependent 
viscoelastic  behavior,  along  with  its  theoretical  basis, 
provides  a  foundation  for  further  uses  of  the  model. 
Demonstrations  by  Bagley,  Torvik  and  others  have  shown  that 
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the  use  of  fractional  derivatives  can  be  used  to  improve 
feedback  control  systems  (6),  conduct  transient  analysis  of 
elastic  and  viscoelastic  configurations  (1),  predict  material 
response  outside  of  experimental  data  (5),  and  various  other 
appl i cat i ons . 
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II.  Fractional  Derivative  Viscoelastic  Model 


The  basis  of  this  research  is  to  model  viscoelastic 
material  with  fractional  derivatives.  In  this  chapter,  the 
use  of  fractional  derivatives  is  developed  and  incorporated 
into  the  complex  modulus  model.  The  effects  of  strain, 
frequency,  and  temperature  are  also  investigated  since  they 
affect  the  viscoelastic  behavior. 

2 . 1  Fract i ona 1  Per i vat i ve  Formulat ion 

In  an  attempt  to  model  viscoelastic  materials  with 
fractional  powers,  the  use  of  the  fractional  derivative  is 
required.  In  a  paper  by  Torvik  and  Bagley  (5:126),  the 
fractional  derivative  is  in  the  form 

t 

.  Da[X(t)] - i -  f  iitzll  dr  ( 

dra  r(l-a)  dt  J  r 

0 


where 

0  <  a  <  1 

r  is  the  gamma  function 

Using  Leibniz's  Rule  (7:282-84)  to  expand  the  integral 
yields : 


Da[x(t) ) 


T( 


_i _ r  x(  t-T )  dT  +  i  x( o j 

1-cO  j  Ta  r ( 1  -a)  ta 


(2) 
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The  goal  of  this  process  is  to  have  a  usable  form  of  the 
fractional  derivative  for  use  in  the  viscoelastic  model 
formulation.  By  working  in  the  Laplace  domain,  this  is 
achieved.  The  Laplace  transformation  is  defined  as 


£ 


oo 

J  e-St  x(t)  dt 
0 


(3) 


Taking  the  Laplace  transformation  of  the  integral  in  equation 
(2)  leads  to 


1  f  x(  t-r  )  dT  _  s/T[ x ( t )  ]  f  e  3T 

T(l-a)  J  t"  r(l-oO  J  t°* 


dr 


(4) 


Using  the  relationship  (8:22) 

n,  _  T(n+1) 

’  =  n+1 


and  letting  n  =  -a,  equation  (4)  reduces  to 


(5) 


T< l-oO 


dr 


=  sa_1  ^s  £  [x(t)]  -  x(0)j  (6) 
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Similarly, 


x(  0 ) 


r(l-a) 


=  s0*"1  x(  0  ) 


(7) 


The  final  result  is 

£  ^Da[x(t)]J-  =  sa  £[x(t>] 


(8) 


The  result  of  equation  (8)  is  that  it  is  possible  to  solve 
problems  when  the  material  is  modeled  by  fractional 
derivatives . 

2 . 2  Complex  Modulus  Model 

For  analysis  of  linear  viscoelastic  behavior,  a  popular 
frequency-dependent  complex  modulus  form  often  used  is 


E*  =  E'  +  IE'  ' 


(9) 


where 

E*  =  complex  modulus 
E'  =  storage  modulus 
E''»  loss  modulus 

The  storage  modulus,  E'(w),  is  the  real  component  of  the 
complex  modulus  and  represents  the  stress  in  phase  with  the 
strain  divided  by  the  strain.  The  storage  modulus  is  a 
measure  of  the  energy  stored  elastically  and  recovered  per 
cycle  in  the  material  (9:41).  The  loss  modulus,  E''(co),  is 
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the  imaginary  component  of  the  complex  modulus  and  represents 
the  stress  90°  out  of  phase  with  the  strain  divided  by  the 
strain.  The  loss  modulus  is  a  measure  of  the  energy  lost  per 
cycle  due  to  viscous  dissipation  (9:42). 

A  common  relationship  between  the  imaginary  and  real 
components  of  the  complex  modulus  is  the  loss  factor,  which 
is  defined  as 

tan  6  =  511  (10) 

E' 

The  loss  factor,  tan  6,  represents  the  energy  loss  of  the 
viscoelastic  material.  This  energy  loss,  described  in  Payne 
and  Scott  (10:18-22),  can  be  used  to  calculate  the  damping 
factor  of  the  material.  To  maximize  the  use  of  the  damping 
characteristics  of  viscoelastic  material,  a  high  loss  factor 
is  desired.  This  occurs  in  the  transition  region  of 
viscoelastic  materials. 

In  order  to  analyze  the  behavior  of  viscoelastic 
material,  the  complex  modulus  needs  to  be  modeled  over  the 
entire  frequency  range.  In  a  series  of  papers  by  Bagley  and 
Torvik  (1),  (3),  (4),  and  (5),  the  frequency-dependent 
complex  modulus  is  represented  in  the  following  form: 

E(  io>)  =  E  +  E  (  io))a  (11) 

o  1 
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Expanding  equation  (11)  yields: 


E(iw)  =  E  +  «aE(i)a 

o  1 


(12) 


The  following  relationships  are  used  to  expand  equation  (12): 


i)«  =  (e1"'2]* 


(13) 


f  in/2")  na  .  .  na 

e  I  =  cos  y~  +  *  sin  2 — 


(  14) 


Substituting  equations  (13)  and  (14)  into  equation  (12) 
results  in  the  following  complex  modulus: 


E(  iu)  =  E  +  <o 

o 


“E‘  [ 


na  .  . 

cos  +  i  si 


na  1 

D5-J 


(15) 


Separating  equation  (15)  into  its  real  and  imaginary 
components  results  in  the  fractional  derivative  model  for  the 
storage  and  loss  moduli.  The  storage  modulus  is  the  real 
component  of  the  complex  modulus  and  is  defined  as 

E'  (w)  =  E  +  wa  E  cos  ^  (16) 

o  1  Z 

Similarly,  the  loss  modulus  is  the  imaginary  component  of  the 
complex  modulus  and  is  defined  as: 

E'  '  (w)  =  ua  E±  sin  ^  (17) 
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The  result  of  equations  (16)  and  (17)  is  that  the  complex 
modulus  for  viscoelastic  material  can  be  modeled  by  the  use 
of  fractional  derivatives.  Since  fractional  derivatives  are 
easier  to  manipulate  in  the  Laplace  domain,  equation  (11)  is 
transformed  to  the  form 


E(s)  =  E  +  E  sa 

O  1 


(18) 


Equation  (18)  is  limited  to  frequencies  below  the  glass 
transition'  region.  While  this  model  predicts  an  exponential 
increase  in  the  modulus  for  increasing  frequencies  as 
expected,  the  model  is  unbounded  and  does  not  level  off  in 
the  glassy  region.  To  extend  this  model  into  the  higher 
frequency  ranges,  another  parameter  is  included  in  the  model. 
The  complex  modulus  now  has  the  form 


E(s) 


E  +  E  sa 

o _ 1 

,  ,  a 

1  +  bs 


(19) 


The  denominator,  1  +  bsa,  causes  the  storage  modulus  to 
plateau  and  the  loss  modulus  to  diminish  at  high  frequencies 
in  the  glassy  region.  The  result  is  that  the  complex  modulus 
can  be  modeled  by  the  use  of  fractional  derivatives  and  the 
selection  of  parameters  E^,  E  ,  b  and  a. 

2 . 3  Viscoelastic  Material  Propert les 

While  equation  (19)  encompasses  a  broad  frequency  range, 
viscoelastic  behavior  is  also  affected  by  strain  and 
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temperature.  The  effects  of  strain,  frequency,  and 
temperature  are  investigated  and  incorporated  into  the 
complex  modulus  model. 

2*3.1  Strain.  Since  the  use  of  the  complex  modulus  is 
limited  to  linear  viscoelastic  behavior,  the  boundary  of  the 
linear  region  is  required.  In  Figure  1,  the  relationship 
between  strain  level  and  modulus  is  shown. 


Figure  1.  Variation  of  Complex  Modulus  with  Strain  Level 

(11:40) 

For  strains  levels  less  than  1  or  2%,  the  complex  modulus 
remains  linear  (12:261).  By  remaining  below  this  strain 
level,  use  of  the  complex  modulus  is  applicable. 
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2.3.2  Temperature . 


In  the  presence  of  significant 


temperature  gradients,  the  variation  of  the  complex  modulus 
with  temperature  must  be  taken  into  account.  The  effect  of 
temperature  variations  is  shown  in  Figure  2. 
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Figure  2.  Variation  of  Complex  Modulus  with  Temperature 

(11:38) 


In  Figure  2,  the  temperature  scale  has  been  divided  into 
three  regions:  glassy,  transition,  and  rubbery.  The 
transition  from  one  modulus  magnitude  to  another  is  an 
inherent  property  of  viscoelastic  material.  While  the 
modulus  is  in  the  glassy  region,  the  physical  temperature  is 
relatively  low.  At  this  low  temperature,  the  material  is 
"frozen"  and  the  only  movement  occurs  by  straining  the 
inter-molecular  bonds.  This  requires  extremely  high  forces 
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and  thus  accounts  for  the  high  modulus.  At  these  low 
temperatures,  the  material  is  somewhat  brittle  and  the 
behavior  resembles  that  of  glass;  hence,  the  glassy  region. 

As  the  temperature  increases,  it  becomes  easier  and  easier 
for  the  molecular  chains  to  slide  over  one  another.  This 
process  occurs  over  a  narrow  temperature  range  in  the 
transition  region  and  is  characterized  by  a  rapid  change  in 
the  coefficient  of  thermal  expansion.  At  higher 
temperatures,  the  modulus  enters  the  rubbery  region  where  an 
equilibrium  is  reached  (12:230-31). 

The  effects  of  temperature  on  the  molecular  chains  is 
also  described  in  Payne  and  Scott  (10:14).  They  describe  a 
"potential  barrier"  that  must  be  overcome  for  the  molecules 
to  slide  over  one  another.  The  probability  of  overcoming 
this  barrier  is  proportional  to  exp(-E/kT),  where  E  is  an 
activation  energy,  k  is  Boltzmann’s  constant,  and  T  is  the 
absolute  temperature.  As  the  temperature  is  reduced,  the 
probability  of  exceeding  the  potential  barrier  is  reduced. 
This  probability  is  exponential  and  corresponds  to  a  rapid 
change  in  the  modulus  with  changing  temperature.  As  the 
temperature  is  further  reduced,  the  probability  of  movement 
is  greatly  reduced  and  the  material  becomes  glass-like.  The 
opposite  effect  occurs  as  the  temperature  is  increased.  The 
end  result  is  that  in  the  glassy  region,  the  modulus  is  high 
and  the  loss  factor  is  low.  In  the  transition  region,  the 
modulus  is  changing  rapidly  and  the  loss  factor  is  high.  The 
rubbery  region  is  characterized  by  a  low  modulus  and  a  loss 
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factor  that  is  lower  than  the  transition  region. 


2.3.3  Frequency .  The  effect  of  frequency  on  the 
complex  modulus  is  quite  similar  to  temperature  effects.  In 
the  glassy  region,  the  frequency  is  high.  As  the  frequency 
is  reduced,  the  modulus  transitions  to  the  rubbery  region. 

The  result  of  changing  frequencies  is  that  increasing 
frequency  increases  modulus  while  decreasing  the  frequency 
decreases  the  modulus. 

2.3.4  Combined  Tempera ture-Frequency  Effect .  In  order 
to  account  for  both  frequency  and  temperature  changes  In  the 
complex  modulus  model,  a  reduced  frequency  principle  is  used. 
Using  the  method  of  reduced  variables  (9:266-73),  the  effects 
of  temperature  variations  can  be  incorporated  into  a 
frequency-dependent  complex  modulus  model.  The  effect  of 
changing  the  temperature  results  in  a  horizontal  shift  of  the 
complex  modulus  along  the  frequency  axis.  Thus,  a  modulus  at 
frequency  <*>  and  temperature  T  would  be  equivalent  to  the 
complex  modulus  at  a  reference  temperature  T^  and  using  a 
shifted  frequency  <*>a  ,  where  a  is  the  shift  factor.  The 

T  T 

complex  modulus  developed  in  section  2.2,  equation  (11), 
would  now  have  the  following  form: 

E(  lu)  =  E  +  E  (  ia>a  )“  (20) 

O  1  T 

The  shift  factor,  aT>  is  temperature  dependent  and  is 
determined  empirically  by  superimposing  modulus  data  at 
various  temperatures  and  frequencies.  The  empirical 
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expression  used  to  incorporate  temperature  effects  was 
developed  by  Williams,  Landel,  and  Ferry,  and  is  generally 
referred  to  as  the  WLF  equation.  The  WLF  equation  provides 
an  expression  for  the  general  curve  of  the  shift  factor.  The 
equation  for  the  temperature  shift  factor  is  described  in 
Ferry  (9:273-80)  and  has  the  form 

c  (T  -  T  ) 

log  a  =  - 2 —  (2i) 

T  c  +  T  -  T 

2  O 


where 

c.=  coefficients  dependent  upon  material 
T  =  temperature 

Tq=  complex  modulus  reference  temperature 

The  WLF  equation  allows  for  the  use  of  complex  moduli  data 
for  temperatures  other  than  the  reference  temperature  of  the 
viscoelastic  material.  The  result  is  that  frequency  and 
temperature  dependency  can  be  incorporated  into  the 
fractional  derivative  complex  modulus  model. 
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Ill . 


Rod  Element  Formulation 


In  this  chapter,  the  formulation  of  the  equations  of 
motion  for  elastic  and  viscoelastic  rods  are  developed  using 
the  complex  modulus  that  is  modeled  with  the  use  of 
fractional  derivatives.  The  solution  technique  incorporates 
the  temperature  dependency  of  the  complex  modulus,  and  a  two 
dimensional  matrix  formulation  is  developed. 

3 . 1  Rod  Element  Formulat ion  for  Constant  Temperature 

In  the  development  of  the  equations  of  motion  for  a 
slender  rod,  a  simple  rod  element  shown  in  Figure  3  is  used. 


+  u2  ( t ) 


Figure  3.  Rod  Element  with  Constant  Properties 

Using  constant  properties  and  a  constant  temperature,  the 
following  differential  equation  of  motion  for  a  slender  rod 
i s  used  (13:45-46): 


EA 


pA 


where 

u  =  longitudinal  displacement 


(22) 
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If  quasi-static  equilibrium  is  assumed,  the  equation  of 
motion  for  the  rod  reduces  to 


*2 

EA  2_y.  =  0  (23) 

0x 

The  result  of  equation  (23)  is  that  the  displacement  field  is 
a  linear  function  of  the  form 

u  ( x  ,  t )  =  u  +  (  u  —  u  )  (  x/L )  (24) 

1  2  1 


The  displacement  field  shown  in  equation 
the  strain  energy  formulation  for  a  rod 
strain  energy  takes  the  form  of 


(24)  is  then  used 
(13:34-37).  The 


in 


U(t) 


(25) 


where 

U(t)  =  strain  energy 


By  using  the  linear  displacement  field  established  in 
equation  (24),  the  strain  energy  for  the  rod  element  is 


L 

U(t)  =  I  (u2  -  2u.u,  +  u2)dx 

J  2L2  2  12  1 

0 


(26) 
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The  final  step  in  the  development  of  a  matrix  form  for  a  rod 
element  is  to  apply  Cast i gl i ano ' s  First  Theorem  to  the  strain 
energy.  Cast i gl i ano ’ s  theorem,  described  in  Saada  (14:453), 
requires  that 


tfU(t) 

0u . 

1 


F. 

l 


(27) 


where 

Fj  =  external  force 


Applying  Cast i gl i ano ' s  theorem  to  the  strain  energy  in 
equation  (26)  results  in  the  following  matrix  formulation  for 
a  slender  rod: 


EA 

L 


(20) 


Equation  (28)  is  the  equation  of  motion  for  elastic  and 
viscoelastic  rods  with  constant  temperature.  The  formulation 
of  the  stiffness  matrix  agrees  with  the  development  in  the 
text  by  Cook  (15:6-9). 

3 . 2  Temperature  Effects 

While  the  displacement  field  in  equation  (26)  is  limited 
to  constant  temperature,  temperature  changes  can  greatly 
affect  the  complex  modulus.  To  increase  the  usefulness  of 
the  viscoelastic  formulation,  the  viscoelastic  stiffness 
matrix  needs  to  be  altered  to  incorporate  temperature 
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variations.  The  strain  energy  for  a  rod  element  used  in 
equation  ( 25 )  is  used  as  a  starting  point  to  construct  a  new 
stiffness  matrix.  The  differential  equation  of  motion  is 
developed  by  applying  Hamilton’s  principle  to  the  strain 
energy  for  a  rod  (neglecting  the  kinetic  energy  due  to  the 
quasi-static  assumption).  This  method  is  described  in 
Meirovitch  (13:42-46).  When  Hamilton's  principle  is  applied 
to  equation  ( 25 ) ,  the  result  is 


t  L 
2 


HAE(s)  ffu'j 

2  UxJ 


^  dx  dt  =  0 


(29) 


t  0 

l 


Applying  the  variation  to  the  integral  yields  the  following 
equation  of  motion: 


AE ( s )  —  =0  (30) 


Young's  modulus,  E(s),  is  a  function  of  frequency  and 
temperature.  For  this  model,  the  temperature  variation  is  a 
function  of  position,  x,  along  the  rod.  Solving  equation 
(30)  for  the  longitudinal  displacement  leads  to 


r  c 

u(x)  =  — - —  dx 

J  AE( s ) 


(31) 


where 


Cj  =  integration  constants 
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For  equation  (31)  to  be  solved,  the  complex  modulus  must  be 
integrable.  In  this  case,  the  temperature  variation  is  a 
function  of  position,  x,  and  this  temperature  change  can  be 
incorporated  by  the  use  of  reduced  variables  described  in 
section  2.3.4.  An  integrable  complex  modulus  that  includes 
the  temperature  dependency  has  the  form 

E(s)  =  E  ( s  )a  (a)a  (32) 

1  T 


Substituting  equation  (32)  into  equation  (31)  results  in 


u  ( x ) 


J 


. _  .  .a.  .  a 

AE  (s)  (a  ) 

1  T 


dx  + 


c 

2 


(33) 


The  modified  complex  modulus  in  equation  (32)  is  limited  to 
use  in  the  transition  region.  While  this  may  appear  to  place 
narrow  constraints  on  the  model,  it  must  be  understood  that 
the  transition  region  is  where  the  damping  capability  of  the 
viscoelastic  material  is  maximized.  Since  the  viscoelastic 
material  is  used  for  its  damping  ability,  it  is  beneficial  to 
remain  in  the  transition  region.  This  transition  region  can 
be  manipulated  by  use  of  frequency  loading  or  temperature 
applications.  Thus,  careful  selection  of  proper  materials, 
temperature  gradients,  and  frequency  loadings  can  shift  the 
transition  region  to  the  operating  region. 

To  integrate  equation  (33),  an  integrable  temperature 
shift  function  needs  to  be  used.  An  exponential  function  is 


21 


used  to  model  the  temperature  shift  factor,  (a^)01.  Examples 
in  Payne  and  Scott  (10:24-33)  and  Jones  (11:40-43)  show  that 
the  temperature  shift  resembles  an  exponential  function.  The 
temperature  shift  is  modeled  as 


(a)"  =  a  exp[ -b( T-T  )] 

T  o 


(34) 


where 

a,b  =  curve  fit  parameters 

=  reference  temperature  for  E(s)  data 
T  =  actual  temperature 


To  determine  the  temperature  distribution,  the  heat  equation 
for  a  rod  is  used  along  with  the  assumption  of  quasi-static 
equilibrium.  The  heat  equation  reduces  to  the  form: 


a2T 

dx2 


0 


(35) 


Equation  (35)  results  in  a  linear  temperature  distribution  as 
a  function  of  position,  x,  along  the  rod.  For  a  basic  rod  in 
Figure  3  with  T(x=0)  =  T  and  T(x=L)  =  Tg,  the  linear 
temperature  function  is 


T( x )  =  T  +  (T-T  ) ( x/L) 

1  2  1 


(36) 
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Substituting  equation  (36)  into  equation  (34)  results  in 


<v 


=  a 


•f 


exp<-b  [  T  -  T  +  (T 


-  Ti)(x/L)]| 


(37) 


This  temperature  shift  function  is  substituted  into  equation 
(33)  and  then  integrated.  The  result  is 


u(x)  =  c4(L)exp  {-*[ T.  -  Tq  +  (T2  -  Ti)(x/L)]J 


[' 


+  I abAE( s ) ( T 


-  v] 


+  c  (38) 


To  solve  for  the  integration  constants,  the  boundary 

conditions  u(x=0)  *  u  and  u(x=L)  =  u  are  used.  Once  the 

1  2 

displacement  field,  u(x),  is  solved,  it  is  then  substituted 
into  the  strain  energy  equation  for  a  viscoelastic  rod.  The 
strain  energy  has  the  form 


U(s) 


AE( s ) ( a  )a 

T 

2 


(  39) 


Substituting  the  results  of  equation  (38)  into  equation  (39) 
and  then  integrating  yields  the  strain  energy  of  the  form: 


U(s)  *  abAE(s)(Tz  -  -  u2>2 


exp  [b(Tt 


(40) 
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The  final  step  is  to  apply  Cast i gl i ano ’ s  theorem  to  the 
strain  energy  in  equation  (40).  The  result  is  the  following 
equation  of  motion  for  a  viscoelastic  rod  that  incorporates 
temperature  effects: 


CAE(s) 


-1 


11  fui(s)'|  pt(s)' 
1J  lu_  (S)J  If  (S). 


(41 ) 


where 


C  =  ab( T  -T 
2 


exp[b(T  -T  )]  -  exp  [b( T  - 


.-V)} 


Equation  (41)  represents  the  viscoelastic  formulation 
incorporating  a  temperature  variation  along  the  rod  into  the 
complex  modulus.  A  special  case  exists  when  the  temperature 
along  the  rod  is  constant  but  different  from  the  reference 
temperature,  Tq,  for  the  complex  modulus  data.  If  a 
constant  temperature  other  than  the  reference  temperature 
exists,  the  temperature  shift  function  in  equation  (37)  is 
reduced  to 


(aT)a  =  a  exp[-b(T  -  To) ]  (42) 

Since  the  temperature  shift  function  in  equation  (42)  is  not 
a  function  of  position,  x,  along  the  rod,  the  displacement 
field  is  linear  in  the  form  of  equation  (24).  The  end  result 
is  that,  for  a  constant  temperature  other  than  the  reference 
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temperature,  the  rod  formulation  of  equation  (28)  is  modified 
by  substituting  the  complex  modulus  of  equation  (32)  and 
using  the  temperature  shift  function  of  equation  (42). 

3 . 3  Matrix  Assembly 

The  final  part  of  the  formulation  involves  the 
incorporation  of  elastic  and  viscoelastic  elements  into  a 
matrix  structure.  Figure  4  shows  the  details  of  a  long 
slender  rod  composed  of  elastic  element  1  and  viscoelastic 
elements  2  and  3. 


-*•  f  j  (  s ) 


-»  f2(s) 


■+  f  3  (  s  ) 


-*■  *A<s) 


-»  Uj  (  s ) 


-*  U2(s) 


->  u3(s) 


■+  U4<S) 


-*  X 


Figure  4.  Rod  of  Elastic  and  Viscoelastic  Elements 

At  this  point,  the  displacement  field  and  the  forces  will  be 
placed  in  the  Laplace  domain.  Work  in  the  Laplace  domain 
will  allow  easier  manipulation  of  the  mathematical 
operations,  especially  in  the  use  of  fractional  derivatives 
that  are  used  to  model  the  stiffness  of  the  viscoelastic 
elements.  Using  the  format  established  in  equation  (28),  the 
matrix  form  for  the  equation  of  motion  for  Figure  4  can  be 
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constructed  as  follows: 


*kl 

0 

0 


-kl 

kl  +  k2 
-k„ 


0 

-k. 


k2  +  k3 


-k. 


0 

0 

■k. 


’  u  x  (  s  ) ' 

’  f^s)' 

u2(  s ) 

f2(s) 

i 

►  =  4 

u^(  s ) 

f  3(s) 

m 

.  V s)- 

.  VS>* 

>  (A3) 


where 


ki  as  stiffness  of  rod  element  i 


At  this  point,  the  stiffness  matrix  is  split  up  into  elastic 
and  viscoelastic  stiffness  matrices  of  the  form 


E(s)  ^Kvj  ju(s)|  +  ^KeJ  ^u(s)^  =  ^f(s)^  ( A4 ) 

The  complex  modulus  has  been  factored  out  of  the  viscoelastic 
stiffness  matrix.  This  formulation  can  be  extended  to 
two-dimensional  form.  The  rod  in  Figure  5  shows  a  general 
conf iguration. 


f  . 


Figure  5.  General  Rod  Configuration 
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Using  the  transformation  matrix  described  in  Cook  (15:25-26), 
the  two  dimensional  rod  formulation  can  be  represented  by  the 
following  form: 


2 

2 

■ 

►  ^ 

C 

cs 

-c 

-cs 

ui 

2 

2 

AE 

cs 

s 

-cs 

-S 

V  . 

X 

2 

2 

1 

L 

-c 

-cs 

c 

CS 

uj 

2 

2 

-cs 

-S 

cs 

s 

V  . 

* 

L  jJ 

(45  ) 


where 

c  =  cos  © 
s  =  sin  © 
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IV.  Example  Problem 


The  purpose  of  this  example  problem  is  to  demonstrate 
the  solution  technique  based  upon  the  use  of  fractional 
derivatives  to  model  viscoelastic  behavior.  The  solution  is 
formulated  with  the  assumption  of  quasi-static  equilibrium. 
This  requires  that  the  operating  frequency  is  low  enough  that 
inertia  effects  are  considered  negligible.  A  typical  space 
structure  in  Figure  6  is  the  basis  for  the  example. 

v . 

.  i 


u 

1 


Figure  6.  Elastic  and  Viscoelastic  Member  Truss 
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The  truss  is  composed  of  elastic  rods  made  of  aluminum  and 
viscoelastic  rods  made  of  neoprene  ( polychloroprene )  rubber. 
The  viscoelastic  rods  are  distinguished  by  their 
cross-hatching.  The  example  will  also  incorporate  a 
temperature  gradient  across  the  truss.  The  node  temperatures 
are  marked  such  that  nodes  1  and  2  operate  at  temperature  TA; 
nodes  3  and  4  at  temperature  Tg;  and  nodes  5  and  6  at 
temperature  T^. 

In  order  to  analyze  the  system  structural  response,  the 
complex  modulus  for  neoprene  rubber  needs  to  be  modeled.  The 
actual  data  for  neoprene  rubber  is  taken  from  a  paper  by 
Madigosky  and  Lee  (16).  One  note  of  caution  is  that  the 
example  problem  uses  approximate  numbers  for  the  complex 
modulus.  The  main  purpose  is  to  demonstrate  the  solution 
technique.  Full  use  of  the  model  can  be  employed  with  actual 
data  for  a  specified  operating  environment.  The  master  curve 
for  neoprene  rubber  that  was  available  only  contained  complex 
modulus  data  in  the  transition  region;  therefore,  this 
problem  will  be  limited  to  that  region.  To  limit  the  complex 
modulus  to  the  transition  region,  the  complex  modulus  is 
modeled  with  the  form  of  equation  (15)  and  setting  the 
parameter  equal  to  zero.  Using  the  format  established  in 
equations  (16)  and  (17)  with  E  equal  to  zero  produces  the 

O 

storage  and  loss  moduli  models  to  estimate  the  fractional 
derivative  complex  modulus  for  neoprene  rubber.  The  master 
curve  referenced  has  data  for  the  storage  modulus  and  the 


29 


loss  factor  (16:348).  To  develop  the  data  for  the  loss 
modulus,  the  relationship  in  equation  (10)  is  used.  With  the 
data  calculated  for  storage  and  loss  moduli,  the  fractional 
derivative  model  parameters  can  be  estimated.  Examination  of 
the  storage  and  loss  modulus  curves  revealed  that  their 
magnitudes  were  nearly  equal  in  the  transition  region.  For 
this  to  occur,  the  arguments  of  sine  and  cosine  in  equations 
(16)  and  (17)  need  to  be  equal.  The  arguments  are  equal  at 
45°;  thus,  the  fractional  power  a  must  be  equal  to  1/2.  The 
final  parameter,  E  ,  is  then  estimated  to  provide  a  best  fit 
fractional  derivative  model.  The  master  curve  for  storage 
modulus  and  the  fractional  derivative  model  is  shown  in 
Figure  7;  and  the  master  curve  for  loss  modulus  and  its 
fractional  derivative  model  is  displayed  in  Figure  8. 

Figures  (7)  and  (8)  were  created  using  the  GRAPHER  software 
program  (17),  and  the  data  points  used  are  contained  in 
Appendix  A.  The  fractional  derivative  model  uses  the 
parameters 

a  =  0.5 

E^  =  1.4  x  10^  N/m2  (sec)0'5 
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YOUNG’S  MODULUS.  LOG  E\  (Pa) 


Figure  7.  Young’s  Modulus,  E' ,  and  Parameter  Curve 
Fit  at  Tq  =  20°C  for  Neoprene  Rubber 
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Figure  8.  Young's  Modulus,  E' ' ,  and  Parameter  Curve 
Fit  at  =  20°C  for  Neoprene  Rubber 
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To  incorporate  temperature  effects  into  the  viscoelastic 
response,  the  complex  modulus  established  in  equation  (32)  is 
used.  The  temperature  shift  function,  (a  )a,  is  calculated 

T 

using  equation  (21).  Using  a  reference  temperature,  T  ,  of 

O 

20°C  and  material  coefficients  c  =  4.09  and  c  =  67,  a  shift 
function  is  determined  using  the  data  tabulated  in  Table  1. 

TABLE  1:  Temperature  Shift  Data  for  Neoprene  Rubber 


T  (  °C  ) 

T  -  T  ( °C) 

G 

/ 

(a  ) 

T 

4.0 

-16.0 

4.38 

7 . 8 

-12.2 

2.85 

12.5 

-7.5 

1.81 

17.0 

-3.0 

1.25 

25 . 8 

5.8 

0.69 

31.0 

11.0 

0.51 

41.8 

21.8 

0.31 

47.0 

27.0 

0.26 

A  temperature  shift  function  is  estimated  by  plotting  (aT)a 
verses  (T  -  Tq).  An  exponential  function  in  the  form  of 
equation  (34)  is  used  to  obtain  a  best  fit.  The  software 
program  GRAPHER  (17)  is  used  to  calculate  a  best  fit 
exponential  curve  and  the  model  and  actual  data  are  plotted 
in  Figure  9.  To  obtain  an  accurate  temperature  shift  model, 
data  points  at  T  -  T  =21.8  and  27°C  were  not  used  in  the 

o 

curve  fitting  routine.  The  reason  for  this  is  that  to  force 
the  curve  through  these  points  created  larger  errors  at  the 
lower  temperature  region,  where  the  shift  factor  becomes 
exponentially  large.  Deleting  these  two  points  allowed  for  a 
more  representative  curve  fit,  and  the  resulting  temperature 
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TEMPERATURE  (T-T0) 

Figure  9.  Temperature  Shift,  (aT>a,  at  *  20*C 

for  Neoprene  Rubber  (Actual  Data:  Table  1) 
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shift  model  still  approximated  the  deleted  points  quite  well. 
Since  the  model  agrees  well  with  the  actual  data,  the 
following  temperature  shift  function  is  used: 

(a  )a  =  1.10144  expC  -0.07839  (T  -  T  )]  (46) 

T  O 

The  complex  modulus  model  Including  frequency  and 
temperature  dependency  is  the  basis  for  analyzing  the  truss 
in  Figure  6.  The  model  assumes  assumes  quasi-static 
equilibrium.  The  complex  modulus  data  for  neoprene  rubber  in 
Figures  (7)  and  (8)  is  in  the  frequency  range  of  100  to 
100,000  Hz  at  a  reference  temperature  of  20°C.  For  the 
inertia  effects  to  be  considered  negligible,  the  transition 
region  needs  to  be  in  the  region  of  1.0  Hz.  This  can  be 
accomplished  by  using  the  temperature  shift  property 
described  in  section  2.3.4.  Since  lowering  the  temperature 
shifts  the  transition  region  to  a  lower  frequency  domain, 
selecting  an  appropriate  temperature  variation  has  the  effect 
of  shifting  the  transition  region  to  a  frequency  range  where 
the  inertia  effects  can  be  considered  negligible.  The  data 
in  Table  2  represents  the  results  of  applying  the  temperature 
shift  function  to  a  range  of  temperatures  at  a  reference 
frequency,  w,  of  1.0  Hz. 
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TABLE  2:  Frequency  Shift  Due  to  Temperature 
Changes  for  Neoprene  Rubber 


T  (-C) 

T  -  T  (°C) 
o 

a 

T 

log  («  aT) 

0 

-20 

27.9 

1.44 

-10 

-30 

133.8 

2.13 

-20 

-40 

641.9 

2.81 

-30 

-50 

3078.6 

3.49 

-35 

-55 

6742.2 

3.83 

-40 

-60 

14765 . 3 

4.17 

-45 

-65 

32336.0 

4.51 

The  result  of  Table  2  is  that  choosing  temperatures  in  the 
range  of  -20  to  -40°C  at  a  frequency  of  1.0  Hz  is  equivalent 
to  operating  at  the  reference  temperature,  T  =  20°C,  in  the 

O 

2  5 

transition  frequency  region  (10  to  10  Hz).  The  lower  limit 
on  the  temperature  range  is  approximately  -50°C,  the  glass 
transition  temperature  for  neoprene  rubber  (18:1564).  The 
truss  will  be  modeled  with  a  temperature  variation 
horizontally  along  the  truss. 

To  determine  the  system  response,  the  equation  of  motion 
developed  in  equation  (44)  is  used  with  the  following  format: 

E<3)  M  {“!=!}  *  [Ke]  {vS’!}  *  {£<3>}  (47> 

To  solve  for  the  system  response,  the  eigenvalues  and 
eigenvectors  need  to  be  determined  from  the  homogeneous 
solution.  The  complex  modulus  format  in  equation  (32)  is 
substituted  into  equation  (47).  Consolidating  the  constants 
into  the  viscoelastic  stiffness  matrix  results  in  the 
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following  homogeneous  equation  of  motion: 


[  ■■  M  *  H 

The  eigenvalues  and  eigenvectors  are  determined  by  setting 
the  determinant  of  the  stiffness  matrix  equal  to  zero.  The 
system  eigenvectors  represent  a  set  of  orthogonal  vectors 
that  span  the  system  response.  Any  general  motion  can  be 
represented  by  a  superposition  of  the  normal  modes  of  the 
system.  Because  of  this  orthogonality  property,  the  system 
eigenvectors  can  be  used  to  uncouple  the  equations  of  motion 
(13:76-81).  This  is  done  by  the  transformation 

{v5s!}  •  c*]  {w<s>}  <49> 


fu(s)\ 
lv(s)J  ' 


{0} 


(48) 


where 

{w(s)>  =  orthogonal  coordinates  of  system 
Substitution  of  this  transformation  into  equation  (47)  yields 


a 

L  1 

s 

M 

Premultiplying  by 

’  •“  M  * 


[*]  *  [*.]  M 


T 

C*3  transforms  equation  (50)  into 

[\s]  ]  {“"»}  •  C*]T  {f(s>} 


(50) 


(51) 
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The  orthogonal  coordinate  vector  (w(s)>  is  normalized  and 
equation  (51)  is  transformed  to 


■*M  *  N\1  ]  [K]'1]  {"<■>}  -  wT  {*<•>}  <52> 


where 


{;<»>}  * 


orthonormal  coordinates  of  system 


Equation  (52)  is  pre-mul t i pi i ed  by  the  matrix 


r  _n1/2 

M 


This  process  results  in  the  diagonal  viscoelastic  stiffness 
matrix  to  reduce  to  the  identity  matrix  while  the  diagonal 
elastic  stiffness  matrix  reduces  to  the  system  eigenvalue 
matrix.  The  reduced  equation  has  the  form 

30  N]  *  [x]  ]  {“<s>}  ■  w  {£<s)}  <53> 

where 


Equation  (53)  is  the  general  form  of  the  equation  of  motion 
for  the  viscoelastic  formulation  developed  in  a  paper  by 
Bagley  and  Calico  (6).  The  notation  is  modified  to  parallel 
the  solution  technique  developed  by  Bagley  and  Calico. 
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Equation  (53)  now  has  the  form 


w(  s ) 


f  (s) 

a  a 
s  +  a 


(54) 


where 

f(s)  =  j>]  (f (s)> 
a  . 

a  =  system  eigenvalues 
Factoring  out  s-a  from  equation  (54)  yields 


w(s)  =  f(s)  s~a 


Using  the  binomial  expansion  (19:347),  equation  (55)  has  the 
form 


To  obtain  the  system  response  in  the  time  domain,  several 
steps  need  to  be  taken.  First,  the  parameter,  s,  is  factored 
out  of  the  expansion.  To  bring  equation  (56)  into  the  time 
domain,  the  convolution  integral  (13:534-35)  is  used.  This 
example  problem  will  use  a  forcing  function  of  a  unit  step 
for  demonstration  purposes.  Using  the  Laplace  transformation 
property  in  equation  (5),  equation  (56)  has  the  following 


1  + 


(I)  ] 


-1 


(55) 
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form  in  the  time  domain: 


t 

w(t)  =  J  D1-a 
0 


X2t2a 

r< 1+a)  +  n l+2a> 


Xt 


(57) 


The  expansion  term  in  equation  (57)  is  the  alpha  order 
Mi ttag-Lef f ler  function.  Equation  (57)  can  be  represented  in 
the  following  reduced  form: 

t 

w(t)  =  J  D1"0*  {  Ea  [  -<*ta)  J  J  pj  dt  (58) 

0 

Equation  (58)  is  integrated  and  has  the  form 

w(t)  =  -X*1  {  Ea  [  -(Xt*)  ]  -  i  }  p]  (59) 

where 

Ea  [  -<““>  ]  -  Y2 

p  =  0 

This  solution  can  then  be  substituted  back  into  the  general 
system  response  equation: 


r<l  +  ap) 


[*]  T  |w ( t )  J 


(60) 
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4 . 1  Example  Setup 


The  system  response  formulation,  equation  (60),  is 
demonstrated  by  comparing  three  variations  to  the  truss  in 
Figure  6.  The  first  example,  case  A,  is  a  truss  composed 
only  of  elastic  members  (aluminum).  This  case  will  serve  as 
a  baseline  to  determine  the  effects  of  temperature  on  the 
damping  characteristics  of  the  viscoelastic  material.  Case  B 
will  use  the  truss  layout  in  Figure  6  with  elastic  (aluminum) 
and  viscoelastic  (neoprene  rubber)  members  at  a  constant 
temperature.  Case  C  will  also  use  the  same  layout  as  case  B 
but  will  have  a  temperature  variation  horizontally  along  the 
truss.  To  solve  for  the  system  response,  the  following 
parameters  are  used: 

A  =  0.001  m2 

L  =  1.0  m  for  horizontal  and  vertical  rods 
F  =  1000  N  v6  (unit  step  forcing  function) 

E  =  7.0  x  104°  N/m2  (18:2211) 

Et  =  1.4  x  106  N/m2  (sec)°'5 
a  -  0.5 
a  «  1.10144 
b  =  0.07839 

The  truss  in  Figure  6  has  six  nodes,  labeled  1-6,  and  has 
nine  degrees  of  freedom.  The  system  response  will  be 
calculated  for  a  unit  step  function  applied  in  the  vertical 
direction  at  node  6,  (Vg),  and  the  response  at  nodes  v3,  v^ , 
Vg ,  and  Vg  will  be  compared  for  all  three  cases. 
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4.1.1  Case  A.  The  response  for  the  elastic  truss  can 


be  found  by  using  the  format  established  in  equation  (28)  and 
inverting  the  stiffness  matrix. 

4.1.2  Case  B.  This  case  was  originally  intended  to  be 
a  constant  temperature  case.  Using  an  operating  temperature 
of  -20*C  at  a  frequency  of  1.0  Hz,  the  complex  modulus  fell 
in  the  beginning  of  the  transition  region  (see  Table  2  and 
Figures  7  and  8).  The  eigenvalues  and  eigenvectors  were 
determined  by  the  use  of  equation  (48).  Of  the  nine 
eigenvalues  calculated,  five  were  repeated  with  a  value  of 
1973.0.  The  existence  of  repeated  eigenvalues  also  resulted 
in  complex  eigenvectors.  To  achieve  real  eigenvectors,  a 
small  temperature  was  imposed  as  follows:  TA  =  -20.00°C; 

T0  =  -20.01°C;  and  Tc  =  -20.02°C.  While  the  temperature 
variation  seems  negligible,  the  slight  change  was  enough  to 
produce  nine  distinct  eigenvalues  and  real  eigenvectors.  The 
eigenvalues  and  eigenvectors  for  case  B  are  included  in 
Appendix  C. 

4.1.3  Case  C.  The  temperature  variation  for  this  case 
is  chosen  to  capture  the  entire  transition  region  of  the 
complex  modulus  for  neoprene  rubber.  Referencing  Table  2, 
the  node  temperatures  were  chosen  as  follows:  T^  =  -20°C; 

Tg  =  -30°C;  and  Tc  =  -40°C.  The  same  procedure  as  case  B  is 
used  to  determine  the  system  response.  The  eigenvalues  and 
eigenvectors  for  case  C  are  included  in  Appendix  D. 


4.2  Results 


To  determine  the  system  response  for  cases  B  and  C,  the 
Mi ttag-Lef f ler  series  needs  to  be  evaluated.  A  FORTRAN 
program  was  created  to  evaluate  the  Mi ttag-Lef f ler  series  for 
each  eigenvalue.  The  program  can  be  found  in  Appendix  E. 
After  the  Mi ttag-Lef f ler  series  is  calculated  for  the  system 
eigenvalues  at  selected  time  increments,  the  system  response 
can  be  determined.  The  matrix  operations  to  determine  the 
system  eigenvalues,  eigenvectors,  and  response  were 
accomplished  by  use  of  the  software  program  MATLAB  (20). 

The  system  response  for  cases  A,  B,  and  C  are  plotted 
for  selected  nodes.  The  vertical  nodes,  v^  (Figure  10), 
v^  (Figure  11),  v 5  (Figure  12),  and  vg  (Figure  13),  were 
chosen  since  their  displacements  were  significantly  affected 
by  the  vertical  unit  step  forcing  function  applied  at  node  6. 
The  time  increments  displayed  in  Figures  10-13  are  defined  as 


follows:  Tq  = 

o 

o 

_9 

sec;  ^  =  10  sec;  T2  = 

CO 

i 

o 

sec ; 

T3  =  10  7  sec; 

T4 

=  10-6  sec;  T5  =  10-5  sec 

'  T 
’  l6 

-4 

=  10  sec; 

and  T?  =  10~3 

sec . 

The  system  response  for 

each 

node  at  each 

time  increment  are  included  in  Appendix  B  for  case  A;  in 
Appendix  C  for  case  B;  and  in  Appendix  D  for  case  C. 

Examination  of  Figures  10-13  reveals  the  effects  of 
temperature  on  the  system  response.  Since  the  complex 
modulus  for  case  C  encompasses  the  entire  transition  region 
due  to  temperature  variations,  its  damping  ability  is 
increased.  This  is  evident  by  noting  the  slower  response 
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CASE  A 
CASE  B 
CASE  C 


RESPONSE  (•4.0406E— 05  m) 


TIME,  T, 

Figure  11.  System  Response  for  Node 


U  5 
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of  case  C.  While  the  response  for  cases  B  and  C  appear  to 
reach  the  static  response  of  case  A,  they  actually  only 
asymptotically  approach  case  A.  The  responses  that  were 
plotted  contained  data  only  to  four  decimal  places.  Actual 
differences  between  the  cases  take  place  past  this  accuracy. 

The  results  of  Figures  10-13  demonstrate  the  ability  to 
determine  the  system  response  of  a  structure  containing  both 
elastic  and  viscoelastic  elements.  This  is  accomplished  by 
incorporating  temperature  and  frequency  dependency  in  the 
complex  modulus  model  with  fractional  derivatives  and 
temperature  shift  functions. 
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V.  Conclusions/Recommendations 


For  nearly  three  hundred  years,  the  fractional 
derivative  has  gone  virtually  unnoticed.  With  the  increased 
emphasis  on  using  viscoelastic  materials  to  capture  their 
damping  and  energy  dissipation  characteristics,  more  accurate 
and  manageable  techniques  are  needed  to  model  the  material 
behavior.  The  fractional  derivative  has  shown  great  promise 
by  its  ability  to  model  the  complex  modulus  data  of  numerous 
viscoelastic  materials  for  up  to  eight  decades  of  frequency 
with  the  determination  of  three  or  four  parameters.  The 
confidence  in  using  the  fractional  derivative  model  has  been 
enhanced  by  the  theoretical  basis  at  the  molecular  level. 

The  cornerstone  of  this  research  was  the  derivation  of  a 
fractional  derivative  model  for  the  complex  modulus  that 
incorporated  temperature  as  well  as  frequency  dependency.  To 
formulate  the  solution,  the  following  requirements  were 
established: 

1.  Formulation  limited  to  longitudinal  motion. 

2.  Quasi-static  equilibrium  was  specifically 
assumed  to  produce  low  frequency  responses. 

3.  Linear  relationships  allowed  the  use  of  Laplace 
transformations . 

4.  Strain  levels  below  1  -  2%  required. 

5.  Temperatures  must  be  below  the  glass  transition 
temperature  for  the  viscoelastic  material. 


The  temperature  and  frequency  dependency  was  built  into  the 
model  due  to  their  inter-relationship.  By  using  the 
principle  of  reduced  variables,  temperature  variations  can  be 
accounted  for  by  applying  a  temperature  shift  function  to  the 
frequency.  The  inclusion  of  temperature  dependency  limited 
the  model  to  the  transition  region  of  the  complex  modulus. 
While  this  removed  the  rubbery  and  glassy  plateau  regions 
from  the  model,  a  very  versatile  model  was  established  for 
the  transition  region.  Since  the  transition  region  is 
characterized  by  a  high  value  for  the  loss  factor, 
maintaining  the  operating  environment  in  this  region 
maximizes  the  damping  ability  of  the  viscoelastic  material. 

The  construction  of  the  temperature  shift  function  by 
use  of  the  WLF  equation  provided  the  necessary  tools  to  solve 
the  example  problem.  This  example  problem  was  chosen  to 
demonstrate  the  ability  to  use  fractional  derivatives  to 
analyze  the  structural  response  of  a  truss  composed  of 
elastic  and  viscoelastic  members  subject  to  a  temperature 
variation  along  the  truss.  The  viscoelastic  rods  were  chosen 
to  be  made  of  a  commercially  available  neoprene  rubber  while 
the  elastic  rods  were  composed  of  aluminum.  In  order  to 
maintain  the  quasi-static  equilibrium  assumption,  the  complex 
modulus  data  needed  to  be  shifted  to  a  lower  frequency  range. 
This  was  accomplished  by  selecting  a  temperature  range  lower 
than  the  reference  temperature  through  the  use  of  the 
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temperature  shift  function.  With  the  appropriate  temperature 
and  materials  selected,  the  fractional  derivative  complex 
model  was  incorporated  into  the  finite  element  formulation  of 
the  equations  of  motion  for  the  truss.  By  working  in  the 
Laplace  domain,  the  solution  was  simplified,  another 
advantage  of  the  fractional  derivative. 

The  example  problem  demonstrated  the  ability  to 
incorporate  the  temperature  and  frequency  dependency  of  the 
viscoelastic  material  for  operating  conditions  in  the 
transition  region.  The  results  of  this  research  and  others 
has  shown  great  potential  for  the  continued  investigation  of 
fractional  derivatives  to  model  viscoelastic  material 
behavior.  The  benefits  such  as  increased  feedback  states, 
simple  but  accurate  models,  robustness,  and  the  ability  to 
predict  material  behavior  outside  of  the  measured  data,  are 
incentives  to  continue  expanding  the  model. 

While  the  fractional  derivative  model  for  the  complex 
modulus  incorporated  temperature  effects  in  the  transition 
region,  future  work  is  needed  to  extend  the  model  to  the 
rubbery  and  glassy  regions.  Follow-on  work  is  also 
suggested  to  evaluate  case  B  at  constant  temperature. 

Solution  techniques  do  exist  for  repeated  eigenvalues, 
although  somewhat  complicated,  and  the  results  in  this 
example  problem  can  serve  as  a  basis  for  comparison. 
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Appendix  A:  Complex  Modulus  and  Model  Data 
for  Neoprene  Rubber  (16) 


Frequency  Actual  Data  Model 

(N/m2)  (N/m2) 


Log  co 

Log  E' 

Log  E* ' 

Log  E' 

Log  E' ' 

2.0 

7.17 

6.96 

7.00 

7.00 

2.1 

7.20 

7.01 

7.05 

7.05 

2.2 

7.23 

7.07 

7.10 

7.10 

2.3 

7 . 27 

7.13 

7.15 

7.15 

2.4 

7.29 

7 .16 

7.20 

7.20 

2.5 

7.32 

7.21 

7.25 

7.25 

2.6 

7 . 37 

7 .27 

7.30 

7.30 

2.7 

7.40 

7 . 31 

7.35 

7.35 

2.8 

7.42 

7.35 

7.40 

7.40 

2.9 

7.47 

7.42 

7.45 

7.45 

3.0 

7.50 

7.46 

7.50 

7.50 

3.1 

7.55 

7.52 

7.55 

7.55 

3.2 

7.59 

7.57 

7.60 

7.60 

3.3 

7.63 

7.62 

7.65 

7 . 65 

3.4 

7.68 

7.68 

7.70 

7 .70 

3.5 

7.72 

7.72 

7.75 

7.75 

3.6 

7.78 

7.79 

7.80 

7.80 

3.7 

7.82 

7.84 

7.85 

7.85 

3.8 

7.88 

7.91 

7.90 

7.90 

3.9 

7.92 

7.96 

7.95 

7.95 

4.0 

7.98 

8.02 

8.00 

8.00 

4 . 1 

8.02 

8.06 

8 . 05 

8.05 

4.2 

8.08 

8.12 

8 .10 

8.10 

4.3 

8.12 

8 .16 

8.15 

8.15 

4.4 

8.20 

8 . 24 

8.20 

8.20 

4.5 

8.28 

8.31 

8.25 

8.25 

4.6 

8.32 

8.35 

8.30 

8.30 

4 . 7 

8.39 

8.41 

8.35 

8 . 35 

4.8 

8.48 

8.49 

8.40 

8.40 

4.9 

8.53 

8.53 

8.45 

8.45 

5.0 

8.58 

8.58 

8.50 

8.50 
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Appendix  B: 


Case  A  Data 


'  v2(t)  ' 

u3(t> 
v  (t) 

«4(t> 
v4(t)  ► 

u5(t) 
v5(t) 

U  -  (  t  ) 

V<>  J 


.0000 
-.3536 
-2.0607 
.  7071 
-1 . 7071 
-.3536 
-4.4749 
1.0607 
-4.4749 
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Appendix  C:  Case  B  Data 


v2(t)  " 
u  3  ( t ) 
v3(t) 

U.  (t) 

v4(t)  . 

u5(t) 
v5(t) 
u6(t) 


13706.4 

0 

0 

0 

0 

0 

0 

0 

0 

0 

4112.0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

284.0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

946.0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1970.6 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1972.5 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1971.9 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1972.7 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1971.2 

.0214 

.2891 

.0439 

.5202 

.1149 

-.8275 

.2877 

.0000 

.0000 

-.1104 

-.0111 

-.2265 

-.0203 

.0656 

-.0153 

.5785 

.5000 

.0000 

.1835 

.7252 

-.4679 

-1.0568 

.0286 

-.8110 

-.2914 

-.5000 

.0000 

.1104 

.0111 

.2265 

.0203 

-.0856 

.0153 

-.5785 

.5000 

.0000 

.2494 

.8751 

-.3327 

-.7863 

.0651 

-.0141 

.5778 

-.5000 

.0000 

-.1548 

.1276 

-.3179 

.2298 

.5396 

.3116 

.0002 

.5000 

1.0000 

.6221 

.5543 

-1.3249 

-.2290 

-.5397 

-.3110 

-.0002 

-.5000 

-1.0000 

.1548 

-.1276 

.3179 

-.2298 

-.5396 

-.3116 

-.0002 

.5000 

1.0000 

.6666 

.4157 

-1.2336 

-.4793 

.6544 

-.5154 

.2879 

-.5000 

-1.0000 

To 

T1 

T2 

T3 

T4 

T  5 

T6 

T7 

.0000 

.0207 

.0550 

.1166 

.1266 

.0479 

.0000 

.0000 

.0000 

-.0362 

-.1018 

-.1022 

-.0837 

.1059 

.3536 

.3536 

.0000 

.1189 

.3354 

.6556 

1.0711 

1.5490 

2.0607 

2.0607 

.0000 

.0126 

.0350 

-.0553 

-.2699 

-.4594 

-.7071 

-.7071 

.0000 

.1355 

.3804 

.6862 

1.0005 

1.3433 

1.7071 

1.7071 

.0000 

-.0557 

-.1586 

-.1739 

-.2380 

.0060 

.3536 

.3536 

.0000 

.2623 

.7528 

1.1644 

1.9082 

3.0262 

4.4749 

4.4749 

.0000 

.0085 

.0251 

-.1413 

-.4691 

-.7131 

-1.0607 

-1.0607 

.0000 

.2818 

.8096 

1.2360 

2.0625 

3.1261 

4.4749 

4.4749 
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'  v2(t)  * 

u  3  (  t  ) 

v3(t) 

vu 

v4(t) 

►  *  1 

Uj(  t> 

v5(t) 

u6(t) 

8011.0 

Appendix  D 

0  0 

0 

Case  C 

0 

Data 

0 

0 

0 

0 

0 

2973.0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

102.0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1107, 

.9  0 

0 

0 

0 

0 

0 

0 

0 

0 

875.4 

0 

0 

0 

0 

0 

0 

0 

0 

0 

414.6 

0 

0 

0 

0 

0 

0 

0 

0 

0 

546.0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1368.0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

624.6 

.0466 

.5073 

-.0164 

-.6225 

-.3639 

-.1763 

.3660 

.0000 

.0000 

-.1259 

.0758 

.2062 

.0970 

-.4319 

.2269 

o 

* 

i 

.5000 

.0000 

.2462 

.7659 

.3665 

.2341 

.9630 

.8962 

-1.3559 

-.5000 

.0000 

.1259 

-.0758 

-.2062 

-.0970 

.4319 

-.2269 

.1540 

.5000 

.0000 

.2944 

.8025 

.2502 

.8693 

.1408 

.6208 

-.9027 

-.5000 

.0000 

-.1455 

.1218 

.3109 

.2615 

.0140 

.2550 

.5298 

.5000 

-1.0000 

.6167 

.5022 

1.2479 

.3077 

.1975 

.1791 

-.9999 

-.5000 

1.0000 

.1455 

-.1218 

-.3109 

-.2615 

-.0140 

-.2550 

-.5298 

.5000 

-1.0000 

.6293 

.4744 

1.1314 

.2325 

.0387 

1.4048 

-.7366 

-.5000 

1.0000 

To 

T1 

T2 

T3 

T4 

T5 

T6 

T7 

.0000 

.0214 

.0521 

.1043 

.1180 

.0301 

.0143 

.0000 

.0000 

-.0284 

-.0660 

-.1266 

-.1743 

-.0247 

.1744 

.3536 

.0000 

.0941 

.2386 

.5312 

.9202 

1.3882 

1.7422 

2.0607 

.0000 

.0118 

.0175 

.0034 

-.1792 

-.3288 

-.5279 

-.7071 

.0000 

.1030 

.2575 

.5598 

.9177 

1.2481 

1.4897 

1.7071 

.0000 

-.0354 

-.0847 

-.2180 

-.2598 

-.2169 

.0834 

.3536 

.0000 

.1703 

.4151 

.9331 

1.3930 

2.1853 

3.3904 

4.4749 

.0000 

.0111 

.0125 

-.0629 

-.2514 

-.4902 

-.7905 

-1.0607 

.0000 

.1749 

.4276 

.965v 

i .45o* 

<!.3v9l 

3.«917 

4.4749 
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Appendix  E:  FORTRAN  Program  to  Calculate 
Mittag-Lef f ler  Series 


PROGRAM  MITLEF 

* 

REAL  LAMDA ,  DENI,  DEN2,  TERM1 ,  TERM2 ,  SUM,  TINC, 

+  MET ( 15  )  ,  TSTART ,  T 

* 

INTEGER  I,  J,  TPTS ,  NMTERM 

★ 

PRINT* ,  'ENTER  NUM  OF  TERMS  IN  M-L  SERIES' 

READ*,  NMTERM 

PRINT*,  'ENTER  STARTING  TIME’ 

READ*,  TSTART 

PRINT*,  'ENTER  TIME  INCREMENT* 

READ*,  TINC 

PRINT*,  'ENTER  NUM  OF  TIME  FOINTS* 

READ*,  TPTS 

PRINT*,  'ENTER  EIGENVALUE* 

READ*,  LAMDA 

★ 

*  ** 

*  **  The  following  loops  read  the  eigenvalue  and  calculate 

*  **  the  Mi ttag-Lef f ler  series.  The  loops  are  designed  to 

*  **  calculate  the  first  two  terms  (TERM1  and  TERM2 )  in 

*  **  the  Mi ttag-Lef f ler  series  and  then  add  further  terms 

*  **  by  multiplying  the  previous  terms  by  appropriate 

*  **  variables.  This  method  determines  the  Mi ttag-Lef f ler 

*  **  series  in  groups  of  two.  The  results  for  each  time 

*  **  increment  are  printed  out  at  the  end  of  the  program 

*  **  in  the  MET  array. 

*  ★* 

★ 

PRINT  20,  LAMDA 

20  FORMAT  (IX,  'LAMDA  =  ',  F8 . 2 ) 

DO  AO  I  =  1 , TPTS 

* 

T  =  TSTART  +  (I-1)*TINC 
T  =  SQRT(T) 

DENI  =  SQRT ( 4 . 0*ATAN (1.0))  *  0.5 
DEN2  =  1.0 

TERM1  =  ( -LAMDA*T ) /DENI 
TERM2  =  ( (LAMDA*T)**2)/DEN2 
SUM  =  1.0  +  TERM1  +  TERM2 
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DO  30  J 


1 , NMTERM 


TERM1  =  TERMl  *  ( ( LAMDA*T ) **2 ) / ( J  +  0.5) 
TERM2  =  TERM2  *  ( ( LAMDA*T ) **2 ) / ( J*2 ) 

SUM  =  SUM  +  TERMl  +  TERM2 
PRINT  28,  SUM 

28  FORMAT  (IX,  *  SUM  =  •,  FI  2. 2) 

30  CONTINUE 

MET ( I )  =  SUM 
PRINT  35,  I, SUM 

35  FORMAT ( IX ,  'TIME  =  ',  13, 2X,  ' POS  *  ',  F12.2) 

40  CONTINUE 

PRINT  55,  (  MET ( I ) ,  I  =  l.TPTS) 

55  FORMAT  ( TPTS ( 1 X , F7 . 4 )  ) 
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